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f^) ' Abstract 

' We discuss a method for computing the generating function for the 

£L|j multiplicity distribution in field theories with strong time dependent ex- 

JL . ternal sources. At leading order, the computation of the generating func- 

' tion reduces to finding a pair of solutions of the classical equations of 

, motion, with non-standard temporal boundary conditions. 

1 Introduction 



X 
b . 

Ever since the advent of high energy colliders, the study of multi-particle final 
states in hadronic collisions has offered the possibility of probing the structure 
of QCD at a deep level. Much attention has been focused on the nature of 
multi-particle production in jets - for a nice review, see Ref. |1I2| . The problem 
is however very general. Theoretical developments in the last couple of decades 
suggest that the bulk of multi-particle production at the highest collider energies 
may be controlled by semi-hard scales. This suggests the possibility that weak 
coupling computations of final states can be compared to experiments. First 
steps in applying this program to understanding high energy nuclear collisions 
at the Relativistic Heavy Collider (RHIC) are very promising We expect 
these considerations to be even more relevant for next generation experiments 
at RHIC and at the Large Hadron Collider (LHC). 

Motivated by these considerations, in a previous paper 0|, we developed 
a general formalism to study the properties of particle production in a field 
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theory coupled to a strong external time-dependent source. We focused there, 
for simplicity, on the case of a (f> 3 scalar field theory. We believe however that 
most of our results are of general validity and can be extended to a gauge theory 
like QCD. 

The relevance of such a field theory coupled to strong sources to the descrip- 
tion of multi-particle final states in high energy hadronic interactions comes 
from the possibility of separating, in a high energy hadron wave-function, par- 
tons with a large longitudinal momentum fraction x from low x partons. This is 
because large x partons are time-dilated and evolve very little over the interac- 
tion time with the other hadron. They are therefore "frozen" color sources for 
partons at smaller x. The latter, in contrast, are dynamical fields over the time 
scales of the interaction process. The Balitsky-Fadin-Kuraev-Lipatov (BFKL) 
evolution equation |5I6| predicts that parton densities grow very rapidly with 
decreasing x. Because this rapid growth leads to large phase-space densities of 
partons in the hadron wave-function, the dynamics of the field that describes 
the small x partons is to a good approximation classical. This separation of the 
degrees of freedom of a high energy hadron into static sources and dynamical 
fields that behave almost classically is the essence of the McLerran-Venugopalan 
model |7I8I9| . The corresponding effective theory is the theory we mimicked in 

In a high energy hadron, the effective theory that results from this separa- 
tion of the degrees of freedom is completely specified by the knowledge of the 
distribution of hard sources. Because physical results should not depend on the 
arbitrary scale which is used to separate the degrees of freedom, the distribution 
of sources must obey a renormalization group equation, known as the JIMWLK 
equation |10lllll2ll3ll4ll5llfill7| . In the limit of a large number of colors and of 
a large nucleus, it takes a much simpler form, the Balitsky-Kovchegov equation 
|18ll9j . This framework is referred to as the Color Glass Condensate (CGC) 
|2()I21I22| . 

In this framework, to compute an observable in a high energy hadronic 
collision, one first finds the appropriate distribution of color sources and sub- 
sequently computes the observable of interest for that particular distribution 
of color sources. Finally, the observable has to be averaged over all possible 
configurations of distributed sources. The former involves solving the evolution 
equation for the distribution of sources. As in 0] we will assume here that the 
distribution of sources is known and focus on the second part of the calculation. 
We will also assume that factorization holds and that any back reaction from 
the produced fields on the sources can be neglected. Finally, we will assume that 
the sources (because of the rapid growth of parton densities with decreasing x) 
are strong sources. This growth is tamed only when the self-interactions of the 
partons become important, a regime which is called "saturation" |23I24I25I26| . 
In this saturation regime, which is natural for a non-Abelian gauge theory at 
sufficiently high energies, the color sources are of the order of the inverse of 
the coupling constant. They are therefore strong sources in a weakly coupled 
theory. 

A crucial feature of the strong source saturation regime is that the calculation 



2 



of any observable in the background provided by such sources is non perturbative 
in the specific sense that one needs to sum an infinite set of Feynman diagrams, 
even at leading order 1 . 

One therefore needs to develop techniques in order to sum the relevant dia- 
grams. In 0] , we showed how the calculation of the average number of produced 
particles can be remapped, at leading order in the coupling constant, into the 
problem of solving the classical equation of motion (EOM) for the field with 
retarded boundary conditions. We also showed that the average multiplicity 
could be computed at next-to-leading order (NLO) by solving in addition the 
equation of motion for a small fluctuation of the field on top of the classical so- 
lution. Again, remarkably, this is an initial value problem with purely retarded 
boundary conditions. The retarded solutions of these two equations of motion 
also allows one in principle to obtain higher moments of the distribution of 
multiplicities. One such quantity is the variance of the distribution of particles. 

This framework was already employed previously to compute gluon produc- 
tion at leading order |37I38I39I4UI41I42| and in the calculation of quark pro- 
duction |43I44| . The latter computation, leading order in quark production, is 
similar to gluon pair production, which contributes to the gluon multiplicity at 
NLO. A recent numerical study of plasma instabilities in heavy ion collisions 
|45I46| is very similar in spirit to the computation of the small fluctuation fields 
necessary to determine NLO contributions to the average multiplicity. One 
should emphasize a crucial aspect of these results: the equations of motion are 
solved with retarded boundary conditions. This allows one to develop straight- 
forward numerical algorithms for Ending these solutions even at next-to-leading 
order. 

In 0] , we introduced a generating function for the distributions of multiplic- 
ities. We only used it there as an intermediate device for finding expressions for 
the moments of the distribution of multiplicities. In fact, determining it would 
allow one to obtain the distribution of probabilities for n-particle final states 
rather than the average multiplicity alone. We shall discuss here the problem of 
computing the generating function for particle production. Expressing the gen- 
erating function as a sum of vacuum- vacuum diagrams in the Schwinger-Keldysh 
formalism |47I48| , we show that the derivative of the generating function is given 
by a sum of diagrams that are very similar to those involved in the calculation 
of the average multiplicity. Further closely examining the set of relevant dia- 
grams at leading order, we show that their sum can be expressed in terms of 
two solutions of the classical equation of motion. However, unlike the average 
multiplicity, the boundary conditions obeyed by these solutions are not retarded 
and instead, one has boundary conditions both at the initial and at the final 
time. This makes the numerical solution of this problem much more difficult. 
We will discuss briefly one possible strategy for its numerical solution. 

The paper is organized as follows. In section |3 we remind the reader of 

1 Note that in collisions involving at least one "dilute" projectile (proton-nucleus collisions 
at moderate energies, or even nucleus-nucleus collisions near the fragmentation region of one 
of the nuclei), the calculation of observables can be carried out analytically 27 28 29 30 31 
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results derived in 0] that are important later in our discussion. In section^ we 
define the generating function and establish a formula for its derivative. This 
formula is valid to all orders in the coupling constant. In section we restrict 
ourselves to leading order; at this order, the derivative of the generating func- 
tion can be expressed in terms of solutions of the classical equation of motion. 
We derive the boundary conditions for these solutions. These are shown to be 
very simple constraints on the fields at times ±00. We briefly discuss a possi- 
ble strategy to compute these solutions numerically. In section we use the 
previous results in order to derive an explicit formula for the variance of the 
number of produced particles, in terms of solutions of the equation of motion 
for small field fluctuations on top of the classical field. Finally, in section [fjl 
we construct a generating function for the distribution of produced energy as 
opposed to the number of particles. Such a generating function is likely more 
relevant to gauge theories that have infrared problems. We show that it can be 
obtained at leading order (in a very similar fashion to the generating function 
for the number) from solutions of the classical EOM. 



2 Reminder of some results of |3j 

In this section, we shall provide a brief synopsis of the results in £0 on computing 
moments of the multiplicity distribution in a field theory with time dependent 
external sources. 



2.1 Model field theory 

The Lagrangian density in our toy model is 

c = &></> - \ m 2 <? - |0 3 + j4> . (1) 

We assume the source j(x) to be time-dependent and of strength 1/g. By that 
we mean that the dimensionless number J d 4 xgj(x) ~ 0(1) when we do the 
power counting for a diagram 2 . 



2.2 Direct calculation of P„ 



When the current j(x) coupled to the fields is time-dependent, non-zero tran- 
sition amplitudes between the vacuum and populated states are allowed. The 
probability P n for the production of n particles is 



|o in ) 



(2) 



It is well known that the transition amplitude (p 1 ■ ■ ■ p n O \it\0in) is the sum of 
all the Feynman diagrams with only sources in the initial state and n particles 

2 Note that a weak coupling approach in this model is only valid for m? > 0. Else, the 
theory has no minimum about which a stable weak coupling expansion can be performed. 



4 




Figure 1: A typical contribution to Pn, the probability to produce a final state 
with 11 particles. The vertical cut delimits the amplitude and its complex 
conjugate. 



in the final state. Note that these Feynman diagrams are built from time- 
ordered propagators. In a field theory coupled to a time-dependent source, one 
should also include the disconnected vacuum- vacuum diagrams in computing the 
probabilities. (This is in sharp contrast to conventional field theories without 
external sources where such diagrams can indeed be ignored because they add 
up to a pure phase in the amplitude and unity in the probabilities.) 

The necessity of keeping track of the vacuum-vacuum diagrams is what 
makes the direct calculation of any of the probabilities P n difficult. In order to 
illustrate the difficulty of the task, we display in figure ^ a typical contribution 
to Pn, the probability of producing 11 particles in the final state. The number 
of disconnected vacuum-vacuum graphs is arbitrary although their sum is the 
exponential of all the connected ones: 



i V = exp 

all 



(3) 



In a field theory in the vacuum, ^ conn ^ 1S a rea ^ number, and exp (i ^ conn V) 
and the similar factor from the complex conjugate amplitude cancel each other. 
However, this is not the case in the presence of a time-dependent source j{x). 
One should also note that it not possible to perform a perturbative expansion 
of the probability P n when the strength of the source is 1/g. Indeed, the lead- 
ing order contribution in this approach corresponds to keeping only the tree 
disconnected factors, but because each of these factors is of order 1/g 2 , any 
truncation is meaningless. 
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2.3 Calculation of the moments 

In 0], we obtained a very compact formula for the probability P n , that reads 



Pn 



iV\j+] p -iV'\j-\ 



] + =J-=J 



(4) 



In this formula, iV[j+] is the sum of all the connected time-ordered vacuum- 
vacuum diagrams, constructed with the source j + , and — iV*[jJ\ is its complex 
conjugate, constructed with the source D\j+,j-] is an operator that acts 
on the sources j+ , j'_ , defined as 



T>\j+,j-} = ~J d 4 xd 4 y G° + _(x,y) (D x + m 2 ){U y + m 2 ) 



Sj+(x) Sj-(y) 



(5) 



where Z is the wave function renormalization constant and G® (x, y) is the 

propagator of the Schwinger-Keldysh |47148| formalism 3 



d 3 p 



(2tt) 3 2E p 



Jp-(x-y) 



(6) 



Eq. Q , although very formal and of no use for calculating P n itself, enables 
one to obtain very simply formulas for the moments. From these, one gets easily 
the average multiplicity 



]+=3-=] 



A crucial observation in 0] was that the product 

e v\j+,j-\ e iV\j+] e -iV[j-] = e iV SK \j+,j-] 



(7) 



(8) 



namely, is the sum of the vacuum-vacuum diagrams of the Schwinger-Keldysh 
formalism, with the source j+ on the upper branch of the time contour, and j_ 
on the lower branch. 

From eq. (Q, one obtains, for (rij 



(n)= d 4 xd 4 y ZG° + _(x,y) It^ (x)T^\y) (x,y) 



J+=J-=J 



(9) 



where and r(+~) are the 1- and 2-point amputated Green's functions in 
the Schwinger-Keldysh formalism: 



rW(*. 

T(+-\x,y) 



D x + fn 2 6iV SK \j+,j-] 
Z 8j± (x) 

D x +m 2 U y + m 2 5 2 iV SK [j+,]-] 



Z 



5j+(x)Sj-(y) 



(10) 



'See appendix A of |1] for a brief reminder on the Schwinger-Keldysh formalism. 
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Diagrammatically, (n) can be represented as 

< n > = CX3 + CD • 

Note that, contrary to the individual probabilities P n , the average multiplicity 
is composed only of connected diagrams. 

A major simplification occurs in this formula at leading order. When j + = 
j— = j, the sum of the 1-point tree diagrams of the Schwinger-Keldysh for- 
malism is equal to the retarded solution of the classical equation of motion, 
with a vanishing boundary condition at xq = — oo. The problem of summing 
all the diagrams that contribute to (n) at leading order therefore maps into 
the problem of finding the solution of a partial differential equation. Moreover, 
because the boundary conditions for this equation are retarded, the problem 
is straightforward to solve numerically. Further, we also showed in 0] that, 
at next-to-leading order, the calculation of (n) only requires one to solve the 
equation of motion for small fluctuations of the field on top of the classical 
solution again, with retarded boundary conditions. This classical solution was 
previously obtained when computing (n) at leading order. 

The simplifications observed for the first moment, the average multiplicity 
(n), are quite generic. For the moment of order p, one can derive a formula 
similar to eq. although these formulas become increasingly cumbersome 

as the order p increases. (The formula for p = 2 was given explicitly in 0], 
and we shall not reproduce it here.) Following the same arguments as those 
for the average multiplicity, all moments can in principle be calculated from 
the retarded solution of the classical EOM plus retarded solutions of the small 
fluctuation EOM. Of course, doing so in practice is prohibitively computing 
intensive with increasing p. 
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3 Generating function: general features 



3.1 Definition 

The generating function is simply defined as 4 

+00 

^■(z)^^P„z", (14) 

where the variable z is possibly complex. 

In 0] , the generating function was purely an intermediate device to simplify 
the derivation of moments of the multiplicity distribution. The generating func- 
tion itself was never computed. However, by calculating T{z), one would have 
access to the complete sequence of probabilities P n ; assuming T(z) is known, it 
is very easy to go back to the probabilities P n by an integration in the complex 
plane 5 : 

P n = -L / T{z) , (15) 

where C is a closed path circling around the origin z = 0. 

Choosing the contour C to be the unit circle, we can obtain P n very effectively 
as a Fourier coefficient 6 of the function T[e %B ) 

Pn = ^~ d9 e- md T(e ld ) . (16) 
£5 Jo 

4 In 2|, we defined it instead as 

+00 

F(x) = Y, Pn e " X ■ ( 12 ) 

n=0 

This alternative definition was chosen because it simplified the calculation of the moments of 
the distribution of multiplicities 



(n p )sY l Pnn p = F( p \0), (13) 

n=0 

to the pth derivative of the generating function. The correspondence between the definition 
in this paper and the one used in |I] is of course F(x) = F(e x ). 

5 From the definition in eq. 1141 . the probabilities P n can also be determined from the 
generating function by taking successive derivatives at z = : P n = J-"' n )(0)/n!. However, 
the interesting values of n are located around the average value (n). Because the typical 
multiplicity is very large in a heavy ion collision, derivatives of very high order would have 
to be calculated for this method to apply. This is difficult to do numerically with a good 
precision. 

B In practice, one should write the Fourier integral as the discrete Riemann sum 

N-l 

p n = lim -J— y e -^W F(e 2 ™W) , 
fc— 

and evaluate the sum by using the "fast Fourier transform" (FFT) algorithm. Note that 
the maximal "frequency" N one uses in the FFT limits the largest multiplicity n which is 
accessible. 
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Figure 2: Distribution of multiplicities for two toy models of the generating 
function !F{z). Solid circles : T\{z) = exp(n(z — 1)). This generates a Poisson 
distribution of average n - we set n = 1000 in the figure. Open triangles : 
T 2 {z) = exp(ri(^p + + + ^y 1 )). In this second model, 25% of the 
particles are produced individually, 25% of the particles come in bunches of 10, 
25% in bunches of 20 and the remaining 25% in bunches of 30. The average n 
is the same as in the first model. 
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The use of this formula to go from the generating function F(z) to the 
distribution of probabilities P n is illustrated in figure with two toy models for 
J-(z). The two models have the same average number of produced particles, 
but differ vastly in how these particles are correlated. As one can see, the effect 
of the correlations in the latter case is to significantly alter the width of the 
multiplicity distribution. It is precisely because the average number (n) alone 
cannot discriminate between these very different production mechanisms that 
one wishes to devise methods to obtain more detailed information about the 
distribution of produced particles. 



3.2 F{z) as a sum of vacuum- vacuum diagrams 

From eqs. (J3J and (|14J) . we obtain simply, 

T{z) = e zV[ ^' ] - ] e lVlj+] e - lV * [3 - ] . (17) 

j+=j-=j 

This expression has a diagrammatic interpretation as the sum of the vacuum- 
vacuum diagrams of the Schwinger-Keldysh formalism, with each off-diagonal 
propagator (H — or — h) weighted by a factor z. This is equivalent to modifying 
the off-diagonal Schwinger-Keldysh propagators by 

r° » y r° 

Lr j r Z U j , 

G°_ + — > z G°_ + . (18) 
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We shall exploit this property later when we propose a method for calculating 
the generating function. For the record, the expressions for all the propagators 
involved in the diagrammatic expansion of the generating function are 



G° ++ (x,y) = 
G°—{x,y) = 
zG° + _(x,y) = z 
zG°_ + (x,y) = z 



d 3 p 



{2tt) 3 2E 1 



d 3 p 



{2tt) 3 2E p 



[o(x° - y °yP< x -y^ + 6{y° - x°)e- ip < x - y A , 



d 3 p 



(2n) 3 2E p 



d 3 p 



(2n) 3 2E p 



a ip-{x-y) 



-ip-(x-y) 



(19) 



where E p = \/p 2 + m 2 is the on-shell energy, and implicitly po — +E p . 
3.3 All-orders result 

The sum of all the vacuum- vacuum diagrams involved in T(z) is the exponential 
of the sum that contains only the connected ones. It can be written as 



e zT>\j + ,j_] e iV\j+] e -iV\j-] = e iV SK [z\j+,j-] 
where V SK [z\j+, j_] is the generalization to z ^ 1 of 

V 3K \j+,j-}=V 3K [l\3+,3-] , 



(20) 



(21) 



encountered earlier in the calculation of (n). The generating function is simply 
obtained from this quantity by equating j + and j_ : 



(22) 



Let us now consider the derivative of the generating function with respect to 
its argument z. We can differentiate directly eq. I|17|) with respect to z, which 
leads to 

(23) 



F(z)=V[j +) j-] ^ SK [z\ j+ ,0-] 



As previously in the calculation of (n), making explicit the action of the operator 
P[j + ,j„], one writes 

T\z) = (? v skW^ J d 4 xd 4 y ZG° + _{x,y)\ i T(+'>{z\x)T(-\z\y) + T( + -*>(z\x,y) 

(24) 



where we have defined 



r (±)( 



z\x) 



D x +m 2 SiV SK [z\j + , j_] 



Z 



Sj± (x) 



T^-\z\x,y) 



D x + m 2 D y + m 2 5 2 iV SK [z\j+, j-} 



Z 



Z 



Sj+(x)Sj-(y) 



(25) 



] + =J-=3 
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Therefore, one can write the derivative as 



Hz) 



d 4 xd 4 y ZG + _(x,y)\T^{z\x)T^-\z\y)+T^ + -\z\x,y)\ . (26) 



This relation is very similar in structure to the formula we derived previously 
in ^ for the average multiplicity (n). In fact, it contains exactly the same 
topologies 7 , 

except that the 1-point and 2-point Greens functions T^(z\x) and \z\x, y) 
must be evaluated with modified - z-dependent - Schwinger-Keldysh rules. It 
is precisely this similarity with (n) which makes the formula very attractive. 
As we shall see shortly, at leading order, T'{z)jT{z) can be expressed in terms 
of solutions of the classical equation of motion. If one succeeds in calculating 
it, T(z) can be subsequently determined from the relation 8 

ln^(z) = ln^(l) + J' dz< . (28) 

The integration constant ln.F(l) can be determined trivially from the unitarity 
condition 

oo 

^(1) = ^P„ = 1. (29) 

n=0 

Thus, 

p{jr*-Q2}. (30) 
4 Generating function at leading order 

An urgent question at this point is whether there exists a practical way to cal- 
culate the ratio T 1 \z) / T \z) , given the strong similarities between the formulas 
for T'{z)jT{z) and for (n). Because there is a well defined algorithm to com- 
pute (n), order by order, in terms of solutions of classical equations of motion 
with retarded boundary conditions, it might be hoped that a similar algorithmic 
procedure exists for T\z)j T{z). 

7 This similarity is not a coincidence. Indeed, 



8 If J 7 ' (z) /J~(z) does not have poles, then we can integrate along any path in the complex 
plane that goes from 1 to z. If on the contrary this quantity has a pole at some zq, then 
In F(z) itself has a logarithmic branch cut starting at zo. The result of the integration now 
depends on the number of windings of the integration path around the pole z$ , an ambiguity 
which amounts to choosing the Riemann sheet in which the point z lies. Note however that 
this does not introduce any ambiguity in T(z) itself. 
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In this section, we explore this issue at leading order. It is important to 
stress at the outset that by "leading order" we refer to the logarithm of the 
generating function, or of its derivative with respect to the variable z. For these 
two quantities, the leading order is the order 1/g 2 . The "next-to-leading order" 
would be the order g°, and so on. Discussing the order of the generating function 
T(z) itself does not make sense because it is evident that the exponentiation 
of the logarithms will mix up terms at all orders. This mixing was previously 
encountered in the probabilities P n themselves 9 . 



4.1 Expression for T'{z)j T{£) at leading order 

At leading order, as for the case of the average multiplicity, only the first term 
in eq. H27J1 contributes. (The second term, having at least one loop, starts at 
order g°.) Moreover, the two Green's functions T^> should be evaluated at 
tree level, and the wave-function renormalization factor Z is equal to one. The 
tree-level 1-point functions 



J4. =J_ =3 



5j±(x) 

are represented diagrammatically by tree diagrams with one external leg: 



*+(z\x) = ]T X - 
+/- 



7 , <m^) = E 

r 



(31) 



(32) 



These graphs are evaluated with Schwingcr-Kcldysh rules wherein the off-dia- 
gonal propagators are multiplied by a factor z (see eq. (|18f) '). 
In terms of these objects, eq. (12611 can be written as 



Hz) 



d 4 xd 4 y G° + _(x, y) {U x + m 2 )(D y + m 2 ) <S> + (z\x)<S>-(z\y) 
d 3 p 



(2tt) 3 2E p 



d 4 x e ip - x (n x + m 2 )<$>+{z\x) 
d A y e~ lp - y (n y + m 2 )<S>^(z\y) 



(33) 



Finally, performing an integration by parts (as in eq. (79) of 0]), we obtain 



T\z) 



Hz) 



d 3 p 



{2it) 3 2E p 



d 3 xe ipx {d° x -iE p )<5>+{z\x) 
d 3 ye- ip -i> (d%+iE p )$_(z\y) 



Xq — + 00 



XQ — — O0 

y =—oo 



(34) 



9 See eq. (39) of where all the numbers b r , as well as the a in the exponential, are of 
order 1 in the leading order approximation 
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Therefore, at tree level, it is sufficient to find the behavior at large times of the 
quantities $±(z|x) to determine T' {z)j T{z). It is important to note here that 
we have to keep the contributions from both large positive and negative times 
in eq. I|34|) . A priori, the fields $±(z|x) are not zero in the two limits. 

Of course, the true difficulty lies in summing all the tree diagrams that 
contribute to eqs. (|32[1 . A possible strategy for performing this resummation is 
to reformulate the problem of summing these diagrams as integral equations. 
In this case, $±(z\x) can be expressed as 10 



d*y{G\ + {x,y)[j{y)- 9 -3>l{z\y) 

-zG° + _(x,y)[j(y)-^_(z\y)]} 
*-(*!*) = iJ*V {zG°_ + {x,y) [j(y) - 9 -& + {z\y) 

(Av)] } • 



-G°__{x,y)[j{y)- 9 -& 



(35) 



Note that in these equations, the off-diagonal propagators G® and G°_ + are 

multiplied by z. In principle, one could solve these equations iteratively, start- 
ing from the 3 = approximation. Note also that $>±(z\x) obey the classical 
equation of motion, 



(□ + m 2 )$±(z\x) + |$|(z|x) = j(x) . 

This is easily seen from eqs. I|35|l and from the identities 

(n x + to 2 ) G° ++ (x, y) = -iS(x - y) , 
(O x + to 2 ) G°__ (x, y) = +i5(x — y) , 
(n x +m 2 )G a _ + (x,y) = (a x +m 2 )G a + _(x,y) 



(36) 



(37) 



4.2 Calculation of $ ± from the classical EOM 

In the previous section, we expressed at leading order (tree level) J r '(z)/J-'(z) in 
terms of a pair of fields, <f>±(z\x). The latter arc specified by the integral equa- 
tions (|35(l . We also pointed out that they are solutions of the classical equation 
of motion. However, for this latter property to be of any use in computing 
<f>±(z\x), we must find the boundary conditions they obey. 

These can be determined by a method reminiscent of the standard proof of 
Green's theorem. Multiply the equation of motion for $ + (z|x) 

(By +TO 2 )<f>+(z|y) = j(y) - ^ 2 + (z\y) , (38) 

10 For z = 1, we have = *_(l|x) = <t>c{x); further, G° ++ - G° + _ = G°_ + - G°__ = 

G° R , where G^ is the free retarded Green's function. The two equations in this case are 
equivalent and identically represent the solution of the classical equation of motion. 
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by G+ + (x, y) on the left and the equation 



G++0, y)(u y +m 2 ) = -i8(x - y) 



(39) 



by <$>+(z\y) on the right. Integrate both the resulting expressions over y, and 
subtract them to obtain 



<f + (z\x) = i I d*yG° ++ (x,y) j(y) - f <tf(*|y) 



+i / d 4 y G° ++ (x,y)[ D v ~ D y ]$+(z\y) 



(40) 



Following a similar procedure for the equation of motion of $_(z|x) and the 
expression in eq. l{3"?|l for G° _(x,y), we obtain 



= i I d 4 yG° + ^x,y) [j(y) - |$ 2 _(z|y) 



+i d 4 y G° + _(x,y)[n y -D y ]<S>-(z\y) 



(41) 



Multiplying this equation by z and subtracting it from eq. (|40p. we obtain 



*+(z\x) 



id 



9^2 

2 



'y[G° ++ {x,y) [j(y) 

-zG° + _(x,y) [j(y)-Z*U*\v)]} 

+i J dSj G° ++ (x,y)[B y - U y ]®+{z\y) 

-i J d 4 y zG° + _(x,y)[D y -Dy]<f-(z\y) . (42) 

The first two lines of this equation correspond to the integral equation (|35|) . 
This means that the sum of the third and fourth lines must vanish. 

We will now show that these terms indeed provide the boundary conditions 
we are looking for. They can be written as 

= Jd 4 y{G° ++ (x,y)[Dy-n y ]^ + (z\y)-zG + _(x,y)[ h v -n y ]$-(z\y)} 

d 4 ydy {G% + (x,y)[ d*- d*]*+{z\y) - zG° + _(x,y)[ d%]^(z\y)} ■ 

(43) 

These terms can be written as an integral over the boundary of space-time. 
Assuming that the fields vanish at infinite distance in the spatial directions 11 , we 

11 Alternatively, as is often the case in numerical computations, one may assume periodic 
boundary conditions in space corresponding to the spatial topology of a torus. The boundary 
integral in the spatial directions can be dropped because a torus has no boundary. 
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are left only with the integration at infinite times, with the temporal boundary 
conditions for $-t(z|ar) at j/o — ioo given by 

jd 3 y {G° ++ (x,y)[d° y ~ 0°]<M*ltf) 

<- -> ~| y°=+oo 

-zG° + _(x,y)[ d° y - d y }$-{z\y)\ =0. (44) 

J y u — — oo 

Note that this relation must be satisfied for any x. 

Proceeding in exactly the same way 12 for §>_(z\x), one obtains 

Jd 3 y {zCP_ + (x,y)[d°y-d°]$ + (z\y) 

-&__(x,v)[d° v - dl]*-(*\v) \ =0- (45) 

J y u — — oo 

The boundary conditions in eqs. (|44|l and l|45[) are greatly simplified by 
expressing $±(;z|x) as a sum of plane waves, 

*-w«)= / (27r faL ( / - +) ( z \ x0 >p) e ~ ip ' w + /- -) ( py p ' x } ■ 

(46) 

In these formulae, the variable po is positive and equal to its on-shcll value 
E p = \/p 2 + m 2 . Because §>±(z\x) does not obey the free Klein-Gordon equa- 
tion, the coefficients functions must themselves depend on time. However, by 
assuming that both the source j(x) and the coupling constant g are switched 
off adiabatically at large negative and positive times, the coefficient functions 
f±^(z\x°,p) become constants in the limit of infinite time. 

Substituting eq. (|46|l in eqs. (|44|l and (|45|l and employing the explicit form 
of the propagators G° e , (x, y) in eqs. I|19|l . these boundary conditions can be 
expressed much more simply as boundary conditions of the f±^\z\x°,p) : 





H* 


x° = 


-oo, p) 


= 0, 








H* 


\x° = 


-oo, p) 


= 0, 








\* 


x° = 


+oo,p) 


= zfi + \z\x° = - 


foo,p) , 




fl- 


\z 


\x° = 


+oo,p) 


= zf { _-\z\x° = - 


-foo,p) . 


(47) 



12 The most notable difference is that the equation obeyed by the propagator Gl_ has the 
opposite sign on its right hand side : 

G°_ _(x,y)(O y +m 2 ) = +i5{x - y) . 
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These relations must be satisfied for each momentum mode p. 

Eq. Q34[l can be rewritten in terms of the coefficient functions fj^\z\xo,p) 
at xq = +00 as 

This equation, the boundary conditions in eqs. (|47|l and the fact that the fields 
$^ obey the classical equation of motion, uniquely determine the generating 
function at leading order. 

4.3 Evaluating eq. (|48|) : practical considerations 

We showed in the previous subsection that the eqs. Q36fl . (|46[) . (|48|) and the 
boundary conditions in eqs. (|47|l completely map the problem of finding the 
generating function at leading order into the problem of finding certain solutions 
of the classical equations of motion. In a sense, this is the analog of what 
was achieved in 0] for the average number of produced particles at leading 
order. The only differences are that the boundary conditions now depend on 
the argument z of the generating function and that the <&± are required to 
satisfy boundary conditions at xq — +00. 

Indeed, for 2 = 1, J-'(l)/J-(l) is nothing but (n). We should therefore re- 
cover the results of 0], namely, that $+(l|a;) and <I>_(l|ir) are equal and that 
they are both the retarded solution of the classical EOM with a vanishing 
initial condition. It is easy to see that this result is implied by the eqs. (|47|) . 
The last two equations tell us that <i>+(l|x) and $_(l|a;) have the same coeffi- 
cient functions at x° = +00, implying that <E> + (l|a;) and $_(l|x) are identical 
everywhere : 

Vx, $+(l|x) = $_(l|x) . (49) 
Likewise, the first two boundary equations tell us that 

<P+(l\x° = -00, x) = *_(l|ai = -oo,x) = . (50) 

Therefore $±(l|a;) is nothing but the retarded solution of the equation of motion 
with a vanishing initial condition. 

Although wc found the boundary conditions for <l>±(z|i£), finding the two so- 
lutions of the classical equation of motion that fulfill these boundary conditions 
is much more complicated than finding the retarded solution. The reason for the 
complication is that the boundary conditions in this case are expressed partly 
at x° — —00 and partly at x° = +00. In practice, this means that the problem 
cannot be solved as an initial value problem starting from some known value at 
.To = —00 and solving the EOM forward in time. Instead, most of the methods 
for solving this kind of problem numerically are "relaxation processes" . Gener- 
ically, one perceives them as algorithms where one introduces some fictitious 
"relaxation time" variable £. The simulation begins at £ = with functions $± 
that satisfy all the boundary conditions (it is easy to construct such fields) but 



g(f) 

Hz) 
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not the equation of motion. These fields are then evolved in £ according to the 
equation 

d 6 $± = (n x + m 2 )$± + ^l-j(x) , (51) 

which admits solutions of the EOM as fixed points. The right hand side of the 
previous equation can in principle be replaced by any function that vanishes 
when $± is a solution of the classical EOM. The freedom to chose this function 
could be used in order to ensure that the fixed point is attractive. Finally, at 
each step in the fictitious time £, we only need to make sure that the updat- 
ing procedure preserves the boundary conditions. A numerical algorithm that 
implements this procedure will be discussed in future. We note that somewhat 
similar techniques have been developed recently to study the non-equilibrium 
real time properties of quantum fields 49 . 



5 Variance at Leading Order 

Although it is numerically challenging to find a pair of solutions of the classical 
EOM that obeys the boundary conditions given in eqs. I|47|) . the results of the 
previous section are nevertheless very useful in deriving an expression for the 
variance of the multiplicity distribution, defined as 



<n 2 ) - (nf 



(52) 



The starting point is to differentiate eq. I|48[) with respect to z, and evaluate the 
result at z = 1. This gives the following relation 



<«>l 



^"(1) - (>'(1)) 
d 3 p 



{2nf2E p 



/"(l|+oo,p)/f''(l| + oo,p) 



(+)'/ 



+/(+)(l|+co,p)/i ^ll + oo.p) 



(53) 



(We used the identity ^"(1) = 1 to get rid of the denominators.) The prime 
symbol / denotes the differentiation with respect to z. Moreover, we simplified 
the formula by taking into account the fact that, at z = 1, the two fields <£>± 
are equal and thus have the same Fourier coefficients; we did not therefore spell 
out their subscripts ±. (Note however that the z-derivatives of these Fourier 
coefficients may differ.) 

The next step is to recognize that the derivatives (with respect to z)- /± ' 
- of the Fourier coefficients of $± are the Fourier coefficients of the derivatives 
<f>' ± . Differentiating the EOM in eq. I|36|l with respect to z, one readily sees that 
these derivatives obey the EOM for small fluctuations on top of the classical 
field, namely, 



(54) 
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where we denote by $(x) = $+(l|a;) = &-(l\x), the retarded solution of the 
classical EOM with null initial conditions. Moreover, one can also differenti- 
ate with respect to z the boundary conditions of eqs. (|47|) to find in turn the 
boundary conditions satisfied by the fields $j_(l|a;). We obtain 





l|z° = 


-oo, p) 


= 0, 










l|x = 


-oo,p) 


= 0, 










l|sc° = 


+00, p) 


= /| +) '( 


l\x° = +00, p) - 


f/ (+) (l|z° = 


foo,p) , 




l\x° = 


+00, p) 


= fL-»( 


l\x° = +00, p) - 


f /(-)(1|.T° = 


-hoo,p) . 



(55) 

Although we still have boundary conditions both at xq = —00 and at xq = +00, 
the problem can be solved as an initial value problem, in contrast, as we have 
seen, with the case of the generating function. This is because the equation of 
motion for small field fluctuations, eq. I|54|l . is linear. The fields $' ± (l|x) that 
obey these boundary conditions can therefore be determined by introducing a 
linear basis for the small fluctuation fields. To do so, first introduce two fields 
77±g(x) obeying the EOM for small fluctuations, that are equal to plane waves 
when xq — > —00, namely, 

D x + m 2 + g$(x) Tj± q {x) = , 

V±q( x ) = e ±lq ' x when x — > —00 . (56) 

This initial condition for T]± q (x) is permitted provided the momentum q is on- 
shell because the classical field $(x) vanishes in the remote past. From the first 
two of eqs. I|55|l . we know that <&+(l|a:) has no positive energy component and 
$'_(l|cc) has no negative energy component at x = —00. Therefore $^(l|a;) 
must be a linear combination of the ?y+q's while $'_(l|x) is a linear combination 
of the r)- q 's : 

The coefficients C± q can be determined from the boundary conditions at xq — 
+00 by decomposing the fields r]± q (x) into positive and negative energy plane 
waves 13 . With the Fourier decomposition 

= / , 2 t P 2E {h { t q \x Q ,p)e-^ x + h { - q \x ,p)^) , 
(x) = J [2 ^ 2Ev {h { ^{xo,p) e-^* + h^(x ,p) e»-»} , (58) 



n- 



13 Although of course, at xo = —00, they contain only one of the two components, they 
acquire the other component during the evolution through the background classical field. 
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we may now rewrite the boundary conditions at xq = +00 in terms of the 
Fourier modes. We obtain, 

/«(l| + oc,p) , 

/<->(l| + oo,p) . 

(59) 

These equations must be satisfied for all the momenta p and are a linear system 
of equations to determine the coefficients C± q . Once they have been found, one 
can express the variance at leading order in terms of these previously introduced 
objects 

° - <»> 1 ,0 = J J^E p Jt$kf q {/<"> (II + oo, p) *S p) C +q 

+fl+\l\ + oo,p)ht](+oo,p)C- q } . (60) 

In conclusion, the variance (at leading order) of the number of produced 
particles, may be computed as follows: 

(i) Find the solution $(a;) of the classical equation of motion, with null initial 
conditions, and calculate its Fourier modes /(±)(+co,p) at large positive 
times. This computation will in principle have already been done when 
evaluating numerically the average multiplicity at leading order. 

(ii) Find the functions r\± q (x) by solving the EOM l|54|) for small field fluctu- 
ations about the classical field $(x), with initial conditions exp(±iq ■ x). 
Usually, one will have discretized the spatial volume, and this step must 
be carried out for each momentum q of the dual lattice. 

(hi) The next step (in analogy with the procedure in (i)) is to calculate the 
Fourier modes h^ q (+oo,p) of the i]± q ' at large positive times. 

(iv) With /(±)(+oo,p) and h± q (+oo,p) in hand, solve the linear system of 
equations (|59|l by inverting a matrix (albeit a large one) to find the coef- 
ficients C± q . 

(v) Finally, one obtains the variance from eq. (|60|) . 

As one can see, the calculation of the variance requires a lot more computational 
work than the calculation of the average multiplicity. However, unlike in the 
case of the generating function itself, all the partial differential equations that 
need to be solved are with retarded boundary conditions. All the other steps are 
"elementary" because they involve only Fourier decompositions or the matrix 
inversion of a linear system of equations. 



19 



d 3 q 
{2n) 3 2E q 



h W q (+oo,p) C- q - h<$ (+00, p) C +q 



/iy (+c»,p) c +q - h y _ q '(+^,p) c- q 



(-)/ 



6 Generating function 

for the energy distribution 

Thus far, we only discussed generating functions for the probabilities for produc- 
ing a given number of particles. However, the distribution of radiated energy 14 
may be better defined in the infrared for gauge theories with massless particles. 

We denote P(E) as the probability density for the energy with P(E)dE 
being the probability of radiating energy between the values E and E + dE. 
One may write P{E) in terms of transition amplitudes as 

P W = E ^ J IT T^Jfe: -E^)|(Pl- Pn out 1 in ) I . (61) 



71=0 



In this formula, Ei is the energy of the particle of momentum p { . 
Replacing the delta function in eq. I|61[) by the identity 



we obtain 



S(E V Ei) = ±- [ d9 e^t^-E) (62) 



-\-oc 

p[E) = k i de ' e ~ ieE ' (63) 



where the generating function T E {&) is defined as 

d 3 Pi 



n— i— 1 



(64) 



It is straightforward to obtain for this generating function a form similar to 
cq. JI2J : 

J~ E {0) = e^fb+J-] e ^b+] e -<vy_] ( (65) 



where we denote 



2?f[j+j_] = iy Arf 4 ./^^!,) (□ x + m 2 )(n„ + m 2 ) . 

(66) 

with 

fife*****- (67) 



14 Although we focus here on the generating function for the distribution of energy, our 
derivation holds for any other quantity that is additive for a set of n particles. 
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This means that the generating function T E (0) is the sum of the vacuum- vacuum 
diagrams in a Schwingcr-Kcldysh formalism in which the off-diagonal propaga- 
tors have been modified as follows : 

< ~ r +- > <J +-,e ' 

G°_ + — G°ll e . (68) 

At leading order, the derivative of lnJ r E (6) can be written in terms of a pair 
of solutions of the classical equation of motion, with the following boundary 
conditions : 

f ( +\d\x° = -oo,p) = 0, 

f { _-\e\x° = -cx>,p) = o, 

f { _ +) (6\x° = +00, p) = e i9E » f { +\6\x° = +oo,p) , 

fi~\e\x° = +00, p) = e ieE " f { _~\e\x° = +00, p) . (69) 

As one can see, calculating the generating function for the energy distribution 
involves solving the classical equation of motion with boundary conditions sim- 
ilar to those for the number distributions. 

7 Conclusions 

In this paper, we discussed the generating function for the distribution of pro- 
duced particles in a field theory coupled to a strong time-dependent source. 
We obtained a general formula for the logarithmic derivative of the generating 
function. At leading order, this formula can be expressed in terms of a pair of 
solutions of the classical EOM. We found that these solutions must obey bound- 
ary conditions both at the initial and final time. Finding numerical solutions 
of a non-linear EOM that obey these boundary conditions is much more diffi- 
cult that solving the EOM with retarded boundary conditions. At present, this 
problem is unsolved; we speculate that numerical "relaxation methods" may be 
applied to solve this problem. 

From the results obtained for the generating function, we also sketched an 
algorithm for calculating the variance of the number of produced particles, which 
only involves partial differential equations with retarded boundary conditions. 

Finally, we showed that the generating function for the distribution of pro- 
duced energy may in principle be calculated by very similar methods. 

In addition to the numerical investigations alluded to in the previous para- 
graphs, several extensions of this work are being considered. One of these is 
to generalize the present study to the situation where one considers particle 
production only in some restricted portion of the phase-space. When particle 
production is forbidden in the complementary part of phase-space, the corre- 
sponding generating function could be used in order to study the precise rela- 
tion of colorless partonic configurations in the color glass condensate framework 
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to experimentally observed rapidity gaps. Another possible extension of this 
work is to derive an evolution equation that would drive the dependence of 
the generating function with the center of mass energy of the collision under 
consideration. 
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